#create data to play with
ID = rep(c("t1","t2","t3","t4","t5"),100)
x = runif(500,min=0,max=100)
y = runif(500,min=0,max=x)
z = rnorm(500,mean=10,sd=3)
tdata = data.frame(ID,x,y,z)

#look at tdata ... see R ref card for more inclusive list
tdata[1:10,] #look at the first 10 rows of data
names(tdata) #gives the column names
summary(tdata) #gives basic statistics of each column
mean(tdata[-1]) #gives the mean of each column
sd(tdata[-1]) #gives the sd of each column
quantile(tdata$y) #get the quantile info of the data
quantile(tdata$y, probs=c(0.025,0.1,0.25,0.5,0.75,0.9,0.975)) #quantiles at specified percentiles
sum(tdata$x) #sums the x column
var(tdata[-1]) #gives the variance/covariance matrix of all but the first column
cor(tdata[-1]) #"pearson" (default), "kendall", or "spearman" correlations
unique(tdata$ID) # get the unique values of ID
levels(factor(tdata$ID)) #same as unique but output as a vector

tmp = tapply(tdata$z, tdata$ID, stem) #stem and leaf plot for each group
sumstat = function(x) {return(c(sum=sum(x), mean=mean(x), var=var(x), sd=sd(x), n=length(x)))}
tmp = tapply(tdata$z, tdata$ID, sumstat) #basic stats for each group

#look at the data to see the relationship between x & y
plot(tdata) #produces matplot
coplot(y~x+z|ID,data=tdata) #conditioning plot
example("coplot") #see different examples of conditioning plots
boxplot(tdata) #boxplot of all data
boxplot(tdata$x ~ tdata$ID)

#Transforming data
tdata$logz = log(tdata$z) #simple log transformation
tdata$sqrtx = sqrt(tdata$x) #square root transformation
tdata[1:10,]

#Advanced transformations - Box Cox example
library(car)
?box.cox
miny = min(tdata$y)
if (miny <=0) {miny = abs(miny)+0.0001} else {miny = 0} # get the min value to ensure everything is > 0
tt = box.cox.powers(na.omit(tdata$y)+miny) #get the start values for tranforming
tdata$boxcoxy = box.cox(tdata$y,start=miny,p=tt$start) #append new column to dataframe
par(mfrow = c(1,2), pty = "s")
hist(tdata$y,freq=F) #plot the original data
curve(dnorm(x, mean=mean(tdata$y), sd=sd(tdata$y)), add=TRUE, lty=2) #add a normal curve
hist(tdata$boxcoxy,freq=F) #plot the transformed data
curve(dnorm(x, mean=mean(tdata$boxcoxy), sd=sd(tdata$boxcoxy)), add=TRUE, lty=2) #add a normal curve
dev.off() #turn off the plot to reset the plotting parameters to default

#linear regression
?lm
lm1 = lm(y~x+z,data=tdata)
lm1 #shows results of lm model
lm1 = lm(y~0+x+z,data=tdata) #fixing the intercept to 0
lm1
str(lm1) #see the attributes associated with the lm
lm1$fitted.values
lm1$coefficients

coef(lm1) #coeficients
effects(lm1) 
residuals(lm1) 
fitted(lm1)
predict(lm1)
confint(lm1) #coefficient confidence intervals
logLik(lm1) #log likelihood of the model
AIC(lm1) #AIC of the model with default k = 2 being -2LL + 2k
AIC(lm1,k=log(nrow(tdata))) #to use BIC or SIC

summary(lm1) #puts out the summary information
summary(summary(lm1)) #defines all the objects associated with summary(lm1)
summary(lm1)$coefficients #get the coefficients and all stats

#glm models
?glm
glm1 = glm(y~x+z,data=tdata)
summary(glm1)

#anova of lm or glm models
?anova
anova(lm1) 

#anova of for group or model comparison
fit1=lm(z~ID,data=tdata)
anova(fit1)
fit2=lm(z~ID+x,data=tdata)
anova(fit1,fit2)

#non linear regressions
?nls
nls1 = nls(y~100*(1-exp(-b*x)),data=tdata,start=list(b=0.1))
plot(tdata$x,tdata$y)
abline(lm(y~x,data = tdata),col = "red",lty=2) #one method of adding a line
lines(1:100,predict(nls1,newdata=list(x=1:100)),col="blue") #second method of adding a line
?SSlogis
